Stereoencephalography (SEEG) is a technique used in drug-resistant epilepsy patients that may be a candidate for surgical resection of the epileptogenic zone. Multiple electrodes are placed using a so-called "frame based" stereotactic approach, in our case using the Leksell frame. In our previous paper "Methodology, outcome, safety and in vivo accuracy in traditional frame-based stereoelectroencephalography" by Van der Loo et al (2017) we reported on SEEG electrode implantation accuracy in a cohort of 71 patients who were operated between September 2008 and April 2016, in whom a total of 902 electrodes were implanted. Data for in vivo application accuracy analysis were available for 866 electrodes.
The goal of the current project is to use a public version of this dataset (without any personal identifiers) to predict electrode implantation accuracy by using and comparing different machine learning approaches.
Pieter Kubben, MD, PhD
neurosurgeon @ Maastricht University Medical Center, The Netherlands
The public dataset contains these variables:
The electrodes are the Microdeep depth electrodes by DIXI Medical.
To the limited extent possible in this case I tried to make these FAIR data and adhere to FAIR guiding principles. In practice this meant I introduced the topic, described my data and created a DOI.
Now let's get started.
In [90]:
# import libraries
import turicreate as tc
import h5py
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
plt.style.use('ggplot')
%matplotlib inline
import warnings; warnings.simplefilter('ignore')
#%xmode plain; # shorter error messages
pd.options.mode.chained_assignment = None
# global setting whether to save figures or not
# will save as 300 dpi PNG - all filenames start with "fig_"
save_figures = False
In [2]:
# load data
electrodes = pd.read_csv('../data/electrodes_public.csv')
In [3]:
# find missing values
nan_rows = sum([True for idx,row in electrodes.iterrows() if any(row.isnull())])
print('Nr of rows with missing values:', nan_rows)
In [4]:
def missing_values_table(df):
mis_val = df.isnull().sum()
mis_val_percent = 100 * df.isnull().sum()/len(df)
mis_val_table = pd.concat([mis_val, mis_val_percent], axis=1)
mis_val_table = mis_val_table.rename(columns = {0 : 'Missing Values', 1 : '% of Total Values'})
return mis_val_table
In [5]:
missing_values_table(electrodes)
Out[5]:
In [6]:
#not used features
electrodes.drop(['EntryX', 'EntryY', 'EntryZ','ScrewLength','PlanningRing','PlanningArc'], axis = 1, inplace = True)
In [7]:
# calculate TPLE and remove entry data from dataframe
electrodes['TPLE'] = np.sqrt(np.square(electrodes['TipX'] - electrodes['PlanningX']) +
np.square(electrodes['TipY'] - electrodes['PlanningY']) +
np.square(electrodes['TipZ'] - electrodes['PlanningZ'])
).round(1)
In [8]:
#plt.figure(figsize=[16, 8])
#sns.distplot(electrodes['TPLE'], hist=False, rug=True)
#sns.distplot(electrodes['TPLE'])
In [9]:
plt.figure(figsize=[16, 8])
sns.kdeplot(electrodes[electrodes['PatientPosition']=='Supine']['TPLE'], label='Supine')
sns.kdeplot(electrodes[electrodes['PatientPosition']=='Prone']['TPLE'], label='Prone', color='y')
Out[9]:
It looks like their distributions come from different samples
In [10]:
print('TPLE Mean Supine: {} (SD {}) \nTPLE Mean Prone: {} (SD {})'.format(
round(electrodes[electrodes['PatientPosition']=='Supine']['TPLE'].mean(),2), round(electrodes[electrodes['PatientPosition']=='Supine']['TPLE'].std(),2),
round(electrodes[electrodes['PatientPosition']=='Prone']['TPLE'].mean(),2), round(electrodes[electrodes['PatientPosition']=='Prone']['TPLE'].std(),2)))
Let's prove if indeed are independet samples with a t-test. If we observe a large p-value, then we cannot reject the null hypothesis of identical means. If the p-value is smaller than alpha, then we reject the null hypothesis of equal means.
In [11]:
alpha = 0.05
t_stat, p = stats.mstats.ttest_ind(electrodes[electrodes['PatientPosition']=='Supine']['TPLE'],
electrodes[electrodes['PatientPosition']=='Prone']['TPLE'])
if p < alpha:
print("p= {} The null hypothesis can be rejected".format(p))
else:
print("p= {} The null hypothesis cannot be rejected".format(p))
So it means the two groups came from diferent populations, and we might trate them different
We do the same with electrode type
In [12]:
print('TPLE Mean Oblique: {} (SD {}) \nTPLE Mean Orthogonal: {} (SD {})'.format(
round(electrodes[electrodes['ElectrodeType']=='Oblique']['TPLE'].mean(),2), round(electrodes[electrodes['ElectrodeType']=='Oblique']['TPLE'].std(),2),
round(electrodes[electrodes['ElectrodeType']=='Orthogonal']['TPLE'].mean(),2), round(electrodes[electrodes['ElectrodeType']=='Orthogonal']['TPLE'].std(),2)))
Computing t-test.
In [13]:
alpha = 0.05
t_stat, p = stats.mstats.ttest_ind(electrodes[electrodes['ElectrodeType']=='Oblique']['TPLE'],
electrodes[electrodes['ElectrodeType']=='Orthogonal']['TPLE'])
if p < alpha:
print("p= {} The null hypothesis can be rejected".format(p))
else:
print("p= {} The null hypothesis cannot be rejected".format(p))
In [14]:
#missing_values_table(electrodes)
In [15]:
#taking into account came from different distributions
df_supine = electrodes[electrodes['PatientPosition']=='Supine']
df_prone = electrodes[electrodes['PatientPosition']=='Prone']
We start with the variable to input DuraTipDistancePlanned, but first take a look to the distribution
In [16]:
dura_sup = df_supine['DuraTipDistancePlanned'].dropna()
dura_pro = df_prone['DuraTipDistancePlanned'].dropna()
In [17]:
# normallity test
def normal_test(vector, alpha = 0.05):
k2, p = stats.mstats.normaltest(vector)
if p < alpha: # null hypothesis: x comes from a normal distribution
print("p = {} The null hypothesis can be rejected".format(p))
else:
print("p = {} The null hypothesis cannot be rejected".format(p))
In [18]:
#electrodes['DuraTipDistancePlanned'].hist()
plt.figure(figsize=[16, 8])
sns.distplot(dura_sup)
sns.distplot(dura_pro)
normal_test(dura_sup)
normal_test(dura_pro)
They are normal distribuited so we might input variables with the statistics of the variable
In [19]:
def input_mean(df, variable):
df_null = df[df[variable].isnull()]
df_not_null = df[df[variable].notnull()]
df_null[variable] = round(df_not_null[variable].mean(),1)
inputed_df = df_null.append(df_not_null).sort_index()
return inputed_df
In [20]:
df_supine = input_mean(df_supine, 'DuraTipDistancePlanned')
In [21]:
df_prone = input_mean(df_prone, 'DuraTipDistancePlanned')
We merge again
In [22]:
electrodes = df_supine.append(df_prone).sort_index()
In [23]:
#missing_values_table(electrodes)
Next variable to input is the number of contacts
In [24]:
#electrodes.insert(0, 'Index', electrodes.index)
In [25]:
contact_null = electrodes[electrodes['Contacts'].isnull()]
contact_not_null = electrodes[electrodes['Contacts'].notnull()]
In [26]:
contact_not_null['Contacts'] = contact_not_null['Contacts'].astype(int)
We need to convert into sframe
In [27]:
sf_contact_not_null = tc.SFrame(data=contact_not_null)
In [28]:
# features for contacts classifier
features = ['PatientPosition', 'ElectrodeType', 'PlanningX', 'PlanningY', 'PlanningZ', 'DuraTipDistancePlanned', 'SkinSkullDistance', 'SkullThickness', 'SkullAngle']
In [29]:
# We train a classifier
train_data, test_data = sf_contact_not_null.random_split(0.85)
clf1 = tc.classifier.create(train_data, target='Contacts', features = features)
In [30]:
clf1
Out[30]:
In [31]:
#clf1.evaluate(test_data)
#clf1.evaluate(train_data)
In [32]:
#now convert the rest of the data
sf_contact_null = tc.SFrame(data=contact_null)
new_contacts = clf1.predict(sf_contact_null[features])
contact_null['Contacts'] = new_contacts
New electrodes table with inputed values
In [33]:
electrodes = contact_null.append(contact_not_null).sort_index()
In [34]:
#missing_values_table(electrodes)
# No missing values anymore
How the categories look like on the two different groups?
In [35]:
contacts_supine = df_supine.groupby('Contacts').count()['TPLE']
contacts_prone = df_prone.groupby('Contacts').count()['TPLE']
In [36]:
plt.figure(figsize=[20, 8])
plt.subplot(1,2,1); plt.xticks(contacts_supine.index)
plt.bar(contacts_supine.index, contacts_supine)
plt.subplot(1,2,2), plt.xticks(contacts_prone.index)
plt.bar(contacts_prone.index, contacts_prone, color='y')
Out[36]:
Now we will remove large outliers (difference between planned and real coord) in the Z-axis as electrode insertion length (depth) is influenced also by other factors (calculations regarding depth which could lead to either too superficial or too deep, but also possible malfixation of the screw cap which may cause loosening of the electrode and hence a more superficial position.. it won't migrate into the depth spontaneously). These are very limited numbers and would too much influence further analysis.
In [37]:
plt.figure(figsize=[16, 8])
coordinates = ['PlanningX','PlanningY','PlanningZ','TipX','TipY','TipZ']
for c in coordinates:
sns.kdeplot(electrodes[str(c)], label=str(c))
Every Axis follows a particular distribution, let's try a normal test
In [38]:
for c in coordinates:
normal_test(electrodes[str(c)])
Not normal distribuited
In [39]:
electrodes.head()
Out[39]:
We need to structure our data properly for further analysis and convert the categorical variables (nominal, ordinal) to the category type.
In [40]:
electrodes.drop('TPLE', axis = 1, inplace = True)
In [41]:
#electrodes.insert(0, 'Index', range(0, len(electrodes), 1))
sf_electrodes = tc.SFrame(data=electrodes)
In [42]:
sf_electrodes.head()
Out[42]:
In [43]:
sf_electrodes.num_rows()
Out[43]:
In [44]:
features_train = ['PatientPosition',
'Contacts',
'ElectrodeType',
'PlanningX',
'PlanningY',
'PlanningZ',
'DuraTipDistancePlanned','SkinSkullDistance',
'SkullThickness',
'SkullAngle']
In [45]:
targets_train = ['TipX', 'TipY', 'TipZ']
In [46]:
train_data, test_data = sf_electrodes.random_split(0.8)
In [47]:
def multi_boosted_tree(dataset, features, targets):
models = []
for target in targets:
reg = tc.boosted_trees_regression.create(dataset, target=target, features = features, verbose=False)
models.append(reg)
return models
In [48]:
MBT = multi_boosted_tree(train_data, features_train, targets_train)
In [49]:
MBT
Out[49]:
In [50]:
MBT[1].get_feature_importance()
Out[50]:
In [51]:
# Evaluate the model and save the results into a dictionary
#results =
for i in range(3):
print(MBT[i-1].evaluate(test_data))
In [52]:
test_data['PredictX'] = [round(i,3) for i in MBT[0].predict(test_data)]
test_data['PredictY'] = [round(i,3) for i in MBT[1].predict(test_data)]
test_data['PredictZ'] = [round(i,3) for i in MBT[2].predict(test_data)]
In [53]:
test_data['TPLE'] = np.sqrt(np.square(test_data['TipX'] - test_data['PlanningX']) +
np.square(test_data['TipY'] - test_data['PlanningY']) +
np.square(test_data['TipZ'] - test_data['PlanningZ'])
).round(1)
test_data['TPLE_Pred'] = np.sqrt(np.square(test_data['PredictX'] - test_data['PlanningX']) +
np.square(test_data['PredictY'] - test_data['PlanningY']) +
np.square(test_data['PredictZ'] - test_data['PlanningZ'])
).round(1)
In [54]:
test_data
Out[54]:
In [91]:
#test_data['TPLE'].show()
In [92]:
#test_data['TPLE_Pred'].show()
In [57]:
rmse_1 = tc.evaluation.rmse(test_data['TPLE'], test_data['TPLE_Pred'])
print('Root Mean Squared Error for Boosted Tree Regression \n{}'.format(rmse_1))
In [58]:
td_supine_tple = test_data[test_data['PatientPosition']=='Supine']['TPLE']
td_prone_tple = test_data[test_data['PatientPosition']=='Prone']['TPLE']
In [59]:
print('TPLE Mean Supine: {} (SD {}) \nTPLE Mean Prone: {} (SD {})'.format(
round(td_supine_tple.mean(),2), round(td_supine_tple.std(),2),
round(td_prone_tple.mean(),2), round(td_prone_tple.std(),2)))
In [60]:
estimators = []
for i in range(len(test_data)):
if test_data['PatientPosition'][i] == 'Supine':
estimators.append(round(td_supine_tple.mean(),2))
else:
estimators.append(round(td_prone_tple.mean(),2))
In [61]:
test_data['TPLE_point'] = estimators
In [62]:
rmse_2 = tc.evaluation.rmse(test_data['TPLE'], test_data['TPLE_point'])
print('Root Mean Squared Error for Point Estimator \n{}'.format(rmse_2))
In [63]:
#test_data.print_rows(172,19)
As a quick glance to deep learning we will apply a Multilayer Perceptron (MLP) for multi-class softmax classification using stochastic gradient descent as an optimizer. The code is borrowed from the official keras Sequential model guide and adapted where needed.
In [64]:
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
In [65]:
#First some preprocessing
X = electrodes[features_train].values
#column 0 and 2 are categorical
y = electrodes[targets_train].values
In [66]:
labelencoder_X_1 = LabelEncoder()
X[:, 0] = labelencoder_X_1.fit_transform(X[:, 0])
labelencoder_X_2 = LabelEncoder()
X[:, 2] = labelencoder_X_2.fit_transform(X[:, 2])
In [67]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state=80)
In [68]:
sc = StandardScaler()
X_trainS = sc.fit_transform(X_train)
X_testS = sc.transform(X_test)
In [89]:
from keras.models import Sequential
from keras.layers import Dense, Activation, Dropout
from keras.optimizers import SGD
from keras.utils.np_utils import to_categorical
from keras.callbacks import ModelCheckpoint
In [70]:
batch_size = 2
epochs = 100
In [71]:
model = Sequential()
model.add(Dense(output_dim = 60, activation = 'relu', input_dim = X_train.shape[1])) #(10,)
model.add(Dense(output_dim = 60, activation = 'relu'))
model.add(Dropout(0.5))
model.add(Dense(output_dim = 3, activation = 'linear')) # dim output (3, ) last layer should have linear activations like no activators
In [72]:
model.summary()
In [73]:
#sgd = SGD(lr=0.01, decay=1e-6, momentum=0.9, nesterov=True)
model.compile(optimizer = 'adam', loss = 'mean_squared_error', metrics = ['accuracy']) #mse
In [74]:
#epoch number, 2 decimals - valaccuracy 2 float
filepath = "../models/weights-nn-{epoch:02d}-{val_acc:.2f}.h5"
checkpoint = ModelCheckpoint(filepath, monitor='val_acc', verbose=0, mode='max')
model.fit(X_trainS, y_train,
batch_size=batch_size, verbose=0,
epochs=epochs,
validation_data=(X_testS, y_test),
callbacks=[checkpoint])
Out[74]:
In [76]:
scores = model.evaluate(X_testS, y_test, batch_size=2)
print('\n\nrmse: {:.2f}\n\n{}: {:.2f}'.format(np.sqrt(scores[0]), model.metrics_names[1], scores[1]))
In [77]:
y_pred = model.predict(X_testS)
We are predicting the real carteasian points, so we want to calculate a "real TPLE" before the operation
In [78]:
#feature_idx, feature_name = [], []
#features = []
#for idx, i in enumerate(electrodes.columns):
# features.append((idx,i))
# feature_idx.append(idx)
# feature_name.append(i)
In [79]:
list_coor = ['PlanningX','PlanningY','PlanningZ','TipX','TipY','TipZ']
X = electrodes[list_coor].values #trying the TPLE calculation
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state=80)
In [80]:
coordenates = pd.DataFrame()
coordenates['PlanningX'] = [i for i in X_test[:,0]]
coordenates['PlanningY'] = [i for i in X_test[:,1]]
coordenates['PlanningZ'] = [i for i in X_test[:,2]]
coordenates['TipX'] = [i for i in X_test[:,3]]
coordenates['TipY'] = [i for i in X_test[:,4]]
coordenates['TipZ'] = [i for i in X_test[:,5]]
coordenates['PredictX'] = [y_pred[i][0] for i in range(len(y_pred))]
coordenates['PredictY'] = [y_pred[i][1] for i in range(len(y_pred))]
coordenates['PredictZ'] = [y_pred[i][2] for i in range(len(y_pred))]
In [81]:
coordenates.head()
Out[81]:
In [82]:
coordenates['TPLE'] = np.sqrt(np.square(coordenates['TipX'] - coordenates['PlanningX']) +
np.square(coordenates['TipY'] - coordenates['PlanningY']) +
np.square(coordenates['TipZ'] - coordenates['PlanningZ'])
).round(1)
coordenates['TPLE_Pred'] = np.sqrt(np.square(coordenates['PredictX'] - coordenates['PlanningX']) +
np.square(coordenates['PredictY'] - coordenates['PlanningY']) +
np.square(coordenates['PredictZ'] - coordenates['PlanningZ'])
).round(1)
In [83]:
coordenates.head(10)
Out[83]:
In [84]:
sf_coordenates = tc.SFrame(data=coordenates)
In [85]:
rmse_3 = tc.evaluation.rmse(sf_coordenates['TPLE'], sf_coordenates['TPLE_Pred'])
print('Root Mean Squared Error for Deep Learning \n{}'.format(rmse_3))
In [86]:
rmseS = [rmse_1, rmse_2, rmse_3]
rsme_accuracy = [round(1 - (i/max(coordenates['TPLE'])),4) for i in rmseS]
rsme_accuracy
Out[86]:
In [88]:
#rmse the smaller the better
y=rsme_accuracy#[round(rmse_1,1),round(rmse_2,1),round(rmse_3,1)];
x=list(range(len(y)))
top = [i+0.1 for i in y]; bot = [i-0.1 for i in y]; inter = list(zip(bot,top))
plt.figure(figsize=(16,8))
plt.errorbar(x,y,yerr=[(top-bot)/2 for top,bot in inter], fmt='o', label='Accuracy', color='red', linewidth=3)
plt.plot(x,y, color='black'); plt.legend(); plt.title('Accuracy of RMSE for TPLE')
plt.text(x[0]+.05, y[0]-.01, 'Multi-outupt Boosted Trees')
plt.text(x[1]+.05, y[1]+.01, 'Group Point Estimate')
plt.text(x[2]-.2, y[2]-.01, 'Deep Learning')
Out[88]:
In [ ]: